# Load data on individual donations
load("./data/r02.fd.bd.all.rdata") # outputs an object called "output" into the environment
donations <- output

# We only want to look at first donation event values from each donor
donors <- donations %>%
    group_by(donor) %>%
    filter(date == min(date)) %>%
    ungroup()

# Load FinDonor demographic data
load("./data/r02ds.donorData.rdata") # outputs an object called "output" into the environment
findonor <- output

# Combine the FinDonor datasets
FinDonor <- left_join(donors, findonor, by = "donor")

# Load THL data
# Sofie: thldalta.rdata contains all five THL cohorts, extract FINRISK97 and Health2000 from the others
load("./data/thldata.rdata")
FinRisk97 <- thldata$fr1997
Health2k <- thldata$h2000

# Remove leftovers
rm(output)
rm(thldata)

## Rename useful stuff
# Ferritin, Self-Reported Health
FinDonor <- rename(FinDonor, SRH = QR17, Menstruation = QR79, Age_float = Age, Age = age)
FinRisk97 <- rename(FinRisk97, Ferritin = FERRITIN, SRH = Q40, Gender = SUKUP, Menstruation = K129, Age = IKA)
Health2k <- rename(Health2k, Ferritin = FERRITIINI, SRH = BA01, Gender = SP2, Menstruation = BD03, Age = IKA2, Menopause = MENOP, APOB = NMR_APOB, APOA1 = NMR_APOA1)

# Make "useful stuff" conform with each other
FinDonor <- FinDonor %>% 
    mutate(SRH = case_when(SRH == "Excellent" ~ 1,
                           SRH == "Very_good" ~ 2,
                           SRH == "Good" ~ 3,
                           SRH == "Satisfactory" ~ 4,
                           SRH == "Poor" ~ 5),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == "regular_period" | Menstruation == "irregular_period") ~ "Women|Pre",
                             Gender == "Women" & Menstruation == "no_period" ~ "Women|Post",
                             TRUE ~ "NA")) # Equates to "else"

FinRisk97 <- FinRisk97 %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Pre",
                             Gender == "Women" & Menstruation == 3 ~ "Women|Post",
                             TRUE ~ "NA"))

Health2k <- Health2k %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Pre",
                             Gender == "Women" & (Menstruation == 3 | Age >= 55) ~ "Women|Post",
                             TRUE ~ "NA"))

# Donation eligibility
# These are both "approximates" in a sense, we don't have all the necessary variables to
# filter thoroughly, and we'll be able to do more filtering on Health2000 than FinRisk97
donor_eligible_h2k <- Health2k %>%
    filter(BMII_PAINO.x >= 50 & BMII_PAINO.x <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 & Age <= 66) %>% # Filter away too young and too old
    filter((B_Hb >= 125 & Gender == "Women") | (B_Hb >= 135 & Gender == "Men")) %>% # Filter by hemoglobin
    filter(BA08 == 0) %>% # filter out people with heart attacks
    filter(BA09 == 0) %>% # filter out people with angina
    filter(BA10 == 0) %>% # cardiac insufficiency / heart failure
    filter(!(BA26 == 1 & ATC_A10A == 1)) %>% # filter out people who are diabetic AND use insulin
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) %>% # rename for ease of use
    mutate(HbA1C = B_GHb_A1C * 10.93 - 23.50)

donor_eligible_fr <- FinRisk97 %>%
    filter(PAINO >= 50 & PAINO <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 & Age <= 66) %>% # Filter away too young and too old
    #filter((HGB >= 125 & Gender == 2) | (HGB >= 135 & Gender == 1)) %>% # DON'T filter by hemoglobin, < 500 values in data
    filter(Q15A != 2) %>% # STEMI, NSTEMI
    filter(Q16A != 2) %>% # Stroke
    # filter(Q38 != 2 & Q38 != 4) %>%  # Insulin treatment (2: just insulin, 4: insulin and a tablet)
    filter(Q17B != 2) %>% # cardiac insufficiency
    filter(Q17C != 2) %>% # angina pectoris
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) # rename for ease of use

# Create useful mastersets
# fer_srh <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "SRH", "Group")], 
#                      Health2k = donor_eligible_h2k[, c("Ferritin", "SRH", "Group")], .id = "Cohort") %>% 
#     mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
#            Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k")),
#            SRH = ordered(SRH, levels = 1:5)) %>%
#     filter(Group != "NA") %>%
#     drop_na()

fer_crp <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "CRP")], 
                         Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "CRP")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    filter(CRP >= 0.01) %>%
    drop_na()
table1 <- as.data.frame(table(fer_crp$Group, fer_crp$Cohort))
table1$CRP <- c(paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1$FER <- c(paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1
##         Var1      Var2 Freq                 CRP                      FER
## 1  Women|Pre FinRisk97 1915  0.77 | (0.4, 1.87)   24.02 | (12.39, 42.53)
## 2 Women|Post FinRisk97  884 1.27 | (0.62, 2.65)   55.79 | (31.06, 93.72)
## 3        Men FinRisk97 2603 0.89 | (0.46, 1.88) 112.06 | (66.46, 182.07)
## 4  Women|Pre  Health2k  943 0.62 | (0.27, 1.81)        28 | (15.2, 48.6)
## 5 Women|Post  Health2k  661  1.02 | (0.38, 2.4)      55.9 | (32.7, 95.5)
## 6        Men  Health2k 1710 0.77 | (0.35, 1.75)   124.7 | (76.67, 193.9)
options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    theme_minimal() + 
    geom_smooth(method = "lm", formula = y ~ x, color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort))

options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(aes(color = Group), alpha = 0.2) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(cols = vars(Cohort)) +
    theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(alpha = 1)))

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
CRP_trld <- 3

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_crp %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_crp %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i],var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))) { # run bootstrap only if needed
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_crp %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_crp %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + guides(linetype = "none")

# FINRISK
# Menstruating Women
frwomenpre5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 5)
frwomenpre15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 15)
frwomenpre30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 30)
frwomenpre50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 50)

# Non-menstruating Women
frwomenpost5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 5)
frwomenpost15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 15)
frwomenpost30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 30)
frwomenpost50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 50)

# Men
frmen5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 5)
frmen15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 15)
frmen30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 30)
frmen50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 50)

# H2K
# Menstruating Women
h2kwomenpre5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 5)
h2kwomenpre15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 15)
h2kwomenpre30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 30)
h2kwomenpre50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 50)

# Non-menstruating Women
h2kwomenpost5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 5)
h2kwomenpost15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 15)
h2kwomenpost30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 30)
h2kwomenpost50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 50)

# Men
h2kmen5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 5)
h2kmen15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 15)
h2kmen30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 30)
h2kmen50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 50)
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") +
    theme(legend.position = "none")

ggplot(data = means_all, aes(x = Ferritin, y = means, group = Group)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group, shape = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(cols = vars(Cohort)) +
    labs(y = "%p") + 
    guides(fill = "none") +
    theme(legend.position = "bottom")

# mastersets for the Supplement
# GlycA
fer_glyca <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "GlycA")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "GlycA")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# HbA1C
fer_hba1c <- bind_rows(Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "HbA1C")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()

# APOB
fer_apob <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOB")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOB")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# APOA1
fer_apoa1 <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOA1")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOA1")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
suptable1 <- as.data.frame(table(fer_glyca$Group, fer_glyca$Cohort))

suptable1$GlycA <- c(paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable1$FER <- c(paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable2 <- as.data.frame(table(fer_hba1c$Group))
suptable2$HbA1C <- c(paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable2$FER <-  c(paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable1
##         Var1      Var2 Freq               GlycA                     FER
## 1  Women|Pre FinRisk97 1977 1.28 | (1.17, 1.42)    23.9 | (12.3, 42.52)
## 2 Women|Post FinRisk97  891 1.38 | (1.25, 1.52)   55.78 | (31.2, 93.21)
## 3        Men FinRisk97 2640 1.38 | (1.25, 1.55)  112.1 | (66.53, 183.4)
## 4  Women|Pre  Health2k 1072  1.1 | (0.96, 1.23)  27.64 | (14.78, 48.25)
## 5 Women|Post  Health2k  704     1.15 | (1, 1.3)      55.6 | (32, 95.25)
## 6        Men  Health2k 1903 1.19 | (1.05, 1.35) 121.8 | (75.72, 189.34)
suptable2
##         Var1 Freq                  HbA1C                    FER
## 1  Women|Pre 1074 31.15 | (28.96, 33.34) 27.62 | (14.72, 48.19)
## 2 Women|Post  705 33.34 | (31.15, 36.61)      55.6 | (32, 95.4)
## 3        Men 1910 34.43 | (32.24, 35.52) 121.8 | (75.8, 189.42)
# TODO: suptable3
# TODO: suptable4
options(scipen = 10000)
ggplot(data = fer_glyca, aes(x = Ferritin, y = GlycA)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "FinRisk97", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_glyca %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_glyca %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "Health2k", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_glyca %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_glyca %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")

options(scipen = 10000)
ggplot(data = fer_hba1c, aes(x = Ferritin, y = HbA1C)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group)) +
    theme(legend.position = "none") +
    labs(y = expression(HbA[1*C]))

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
HbA1C_trld <- 42

if (!file.exists(paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_hba1c %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_hba1c %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_hba1c %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))}

means_all <- means_health2k
means_all$Cohort <- c(rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group)) +
    labs(y = "%p") + theme(legend.position = "none")

options(scipen = 10000)
ggplot(data = fer_apob, aes(x = Ferritin, y = APOB)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

options(scipen = 10000)
ggplot(data = fer_apoa1, aes(x = Ferritin, y = APOA1)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
APOB_trld_men <- 1.5 # for males [0.6, 1.5]
APOB_trld_women <- 1.3 # for females [0.6, 1.3]

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_apob %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_apob %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_apob %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_apob %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
APOA1_trld_men <- 1.1 # for males [1.1, 2.0]
APOA1_trld_women <- 1.2 # for females [1.2, 2.3]

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_apoa1 %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_apoa1 %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_apoa1 %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_apoa1 %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))
ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")